Vortex pairing in two-dimensional Bose gases 
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Recent experiments on ultracold Bose gases in two dimensions have provided evidence for the 
existence of the Berezinskii-Kosterhtz-Thouless (BKT) phase via analysis of the interference between 
two independent systems. In this work we study the two-dimensional quantum degenerate Bose gas 
at finite temperature using the projected Gross-Pitaevskii equation classical field method. While 
this describes the highly occupied modes of the gas below a momentum cutoff, we have developed 
a method to incorporate the higher momentum states in our model. We concentrate on finite- 
sized homogeneous systems in order to simplify the analysis of the vortex pairing. We determine 
the dependence of the condensate fraction on temperature and compare this to the calculated 
superfiuid fraction. By measuring the first order correlation function we determine the boundary of 
the Bose-Einstein condensate and BKT phases, and find it is consistent with the superfiuid fraction 
decreasing to zero. We reveal the characteristic unbinding of vortex pairs above the BKT transition 
via a coarse-graining procedure. Finally, we model the procedure used in experiments to infer system 
correlations [Hadzibabic et al.. Nature 441, 1118 (2006)], and quantify its level of agreement with 
directly calculated in situ correlation functions. 

PACS numbers: 03.75.Hh, 03.75.Lm 



I. INTRODUCTION 

At low temperatures a three-dimensional (3D) Bose 
gas can undergo a phase transition to a Bose-Einstein 
condensate. In contrast, thermal fluctuations prevent 
a two-dimensional (2D) Bose gas from making a phase 
transition to an ordered state, in accordance with the 
Mermin-Wagner-Hohenberg theorem [1, 2]. However, the 
2D Bose gas supports topological defects in the form 
of vortices, and in the presence of interactions can in- 
stead undergo a Berezinskii-Kosterlitz-Thouless (BKT) 
[3-5] transition to a quasi-coherent superfiuid state. The 
BKT transition was first observed in liquid helium thin 
films [6], however, more recently, evidence for this transi- 
tion has been found in dilute Bose gases [7-11] (also see 
[12, 13]). 

Ultracold gases have proven to be beautiful systems 
for making direct comparisons between experiment and 
ah initio theory. Experiments in the 2D regime present a 
new challenge for theory as strong fluctuations invalidate 
mean-field theories (e.g. see [5, 14-20]), and only recently 
have quantum Monte Carlo [21, 22] and classical field (c- 
field) [23-25] methods been developed that are directly 
applicable to the experimental regime. 

In this paper we study a uniform Bose gas of finite 
spatial extent and parameters corresponding to current 
experiments. To analyze this system we use the pro- 
jected Gross-Pitaevskii equation (PGPE), a c-field tech- 
nique suited to studying finite temperature Bose fields 
with many highly occupied modes. We develop a tech- 
nique for extracting the superfiuid density based on linear 
response properties, and use this to understand the re- 
lationship between superfiuidity and condensation in the 
finite system. 



With this formalism we then examine two important 
applications: First, we provide a quantitative validation 
of the interference technique used in the ENS experiment 
to determine the nature of two-point correlation in the 
system. To do this we simulate the interference pattern 
generated by allowing two independent 2D systems to 
expand and interfere. Then applying the experimental 
fitting procedure to analyze the interference pattern we 
can extract the inferred two point correlations, which we 
can then compare against the in situ correlations that 
we calculate directly. Second, we examine the correla- 
tions between vortices and antivortices in the system to 
directly quantify the emergence of vortex-antivortex pair- 
ing in the low temperature phase. A similar study was 
made by Giorgetti et al. using a semiclassical field tech- 
nique [26] . We find results for vortex number and vortex 
pair distributions consistent with their results, and we 
show how a coarse graining procedure can be used to 
reveal the unpaired vortices in the system. 

We now briefiy outline the structure of this paper: In 
Sec. II we review the 2D Bose gas and relevant BKT 
physics. In Sec. Ill we outline the c-field technique and 
how it is specialized to describing a uniform, but finite, 
2D Bose gas. In Sec. IV we present our main results, 
before concluding. 



II. FORMALISM 

Here we consider a dilute 2D Bose gas described by the 
Hamiltonian 



Xl/j^'f/'^-i/'V'j (1) 
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where m is the atomic mass, x = {x,y), and ijj = V'(x) is 
the bosonic field operator. 

We take the two-dimensional geometry to be realized 
by tight confinement in the z direction that restricts 
atomic occupation to the lowest z mode. The dimen- 
sionless 2D coupling constant is 



corrections 



9 = 



(2) 



with ttz the spatial extent of the z mode ^ and a the s- 
wave scattering length. We will assume that Uz ^ a so 
that the scattering is approximately three-dimensional 
[27], a condition well-satisfied in the ENS and NIST ex- 
periments [7-9, 11]. For reference, the ENS experiment 
reported in [8] had g 0.15, whereas in the NIST exper- 
iments g w 0.02 [11]. 

In contrast to experiments we focus here on the uni- 
form case; no trapping potential in the xy plane is con- 
sidered. We perform finite sized calculations correspond- 
ing to a square system of size L with periodic boundary 
conditions. Working in the finite size regime simplifies 
the simulations and is more representative of current ex- 
periments. We note that the thermodynamic limit cor- 
responds to taking L — >■ oo while keeping the density, 
n = (V'^V')i constant. 



A. Review of BKT physics 

The BKT superfluid phase has several distinctive char- 
acteristics, which we briefly review. 



1. First order correlations 

Below the BKT transition the first-order correlations 
decay according to an inverse power law: 



9 



(i; 



(x,x') 



(3) 



Systems displaying such algebraic decay are said to ex- 
hibit quasi-long-range order [28]. This is in contrast to 
both the high temperature (disordered phase) in which 
the correlations decay exponentially, and long-range or- 
dered case of the 3D Bose gas in which g^^^ const, for 
llx — x'll — >■ oo. 



2. Superfluid density 

Nelson and Kosterlitz [29] found that the exponent of 

the algebraic decay is related to the ratio of the super- 
fluid density and temperature. To within logarithmic 



a{T) 



(4) 



where ps is the superfluid density and A = h/ \/2iTmkBT 
is the thermal de Broglie wavelength. Furthermore, 
Nelson and Kosterlitz showed that this ratio converges 
to a universal constant as the transition temperature, 
Tkt, is approached from below: limy_^y- a{T) = 1/4 

(i.e., Ps}? = 4). Thus, the superfluid fraction undergoes 

a universal jump from ps{T^r^) = to ps(Tj^rp) = 4/A^ 
as the temperature decreases through Tkt- 



3. Vortex binding transition 

Another important indicator of the BKT transition is 
the behavior of topological excitations, which are quan- 
tized vortices and antivortices in the case of a Bose gas. A 
single vortex has energy which scales with the logarithm 
of the system size. At low temperatures this means that 
the free energy for a single vortex is infinite (in the ther- 
modynamic limit), and vortices cannot exist in isolation. 
As originally argued in [4], the entropic contribution to 
the free energy also scales logarithmically with the sys- 
tem size, and will dominate the free energy at high tem- 
peratures allowing unbound vortices to proliferate. This 
argument provides a simple estimate for the BKT tran- 
sition temperature. 

Although unbound vortices are thermodynamically un- 
favored at T < Tkt, bound pairs of counter-rotating 
vortices may exist since the total energy of such a pair 
is finite ^. This leads to a distinctive qualitative char- 
acterization of the BKT transition: as the temperature 
increases through Tkt pairs of vortices unbind. 



4- Location of the BKT transition in the dilute Bose gas 

While the relation ps(Tj^rj,) = 4/A^ between the super- 
fluid density and temperature at the transition is univer- 
sal, the total density, n, at the transition is not. General 
arguments [30-32] suggest that the transition point for 
the dilute uniform 2D Bose gas is given by 



(nA^)KT = In 



(5) 



where ^ is a constant. Prokofev, Ruebenacker and Svis- 
tunov [14, 15] studied the homogeneous Bose gas us- 
ing Monte Carlo simulations of an equivalent classical 
(j)* model on a lattice. Using an extrapolation to the 



^ For example, for tight harmonic confinement of frequency U)z we 
have az = ^/filrnxoz- 



^ The vortcx-antivortex pair energy depend on the pair size rather 
than the system size. 
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infinite-sized system, they computed a value for the di- 
mensionless constant, ^ = 380 ± 3. By inverting Eq. (5), 
we obtain the BKT critical temperature for the infinite 
system 



KT 



mkB In {^fi^/mg) ' 



(6) 



We use the superscript oo to indicate that this result 
holds in the thermodynamic limit. 



III. METHOD 

A. c-fleld and incoherent regions 

We briefly outline the PGPE formalism, which is de- 
veloped in detail in Rcf. [33]. The Bose field operator is 
split into two parts according to 



■^(x) = Vc(x) +'!/'i(x), 



(7) 



where ipc is the coherent region c-field and tpj is the inco- 
herent field operator (sec [33]). These fields are defined as 
the low and high energy projections of the full quantum 
field operator, separated by the cutoff wave vector K. 
In our theory this cutoff is implemented in terms of the 
plane wave eigenstates {</Jn(x)} of the time-independent 
single particle Hamiltonian, that is, 



k„ = -n, 



(8) 
(9) 



with n = {hx, riy) ^1?. The fields are thus defined by 

V'c(x) = ^ c„(p„(x), (10) 
nec 



where the (in are Bose annihilation operators, the Cn are 
complex amplitudes, and the sets of quantum numbers 
defining the regions are 



C = {n: ||k„|| </^}, 
I = {n: ||k„|| >/^}. 

1. Choice of C region 



(12) 
(13) 



In general, the applicability of the PGPE approach to 
describing the finite temperature gas relies on an appro- 
priate choice for K, so that the modes at the cutoff have 
an average occupation of order unity. In this work we 
choose an average of five or more atoms per mode using 
a procedure discussed in appendix A. This choice means 



that all the modes in C are appreciably occupied, justi- 
fying the classical field replacement a^^ ^ c-^- In contrast 
the I region contains many sparsely occupied modes that 
are particle-like and would be poorly described using a 
classical field approximation. Because our 2D system is 
critical over a wide temperature range, additional care is 
needed in choosing C. Typically strong fiuctuations oc- 
cur in the infrared modes up to the energy scale h^gn/m. 
Above this energy scale the modes are well described by 
mean- field theory (e.g. see the discussion in [14, 34]). For 
the results we present here, we have 



2m 



h^g 
—1 

m 



(14) 



for simulations around the transition region and at high 
temperature. At temperatures well below Tkt, the re- 
quirement of large modal occupation near the cutoff com- 
petes with this condition and we favor the former at the 
expense of violating Eq. (14). 



2. PGPE treatment of C region 



The equation of motion for tpc is the PGPE 



ti^Vi , ti^g^ r, , ,9 , -1 



where the projection operator 

Vc{F{^)} ^ E / rf'xV:(x')F(x'), 



(15) 



(16) 



nGC 



formalizes our basis set restriction of tpc to the C region. 
The main approximation used to arrive at the PGPE is 
to neglect dynamical couplings to the incoherent region 
[35]. 

We assume that Eq. (15) is ergodic [36], so that the 
microstatcs {il>c \ generated through time evolution form 
an unbiased sample of the equilibrium microstates. Time 
averaging can then be used to obtain macroscopic equi- 
librium properties. We generate the time evolution by 
solving the PGPE with three adjustable parameters: (i) 
the cutoff wave vector, K, that defines the division be- 
tween C and I, and hence the number of modes in the C 
region; (ii) the number of C region atoms, A^c; (fii) the 
total energy of the C region, Ec ■ The last two quantities, 
defined as 



Ec 



-I 
I 



-^ + ^l^c|jVc, (17) 

(18) 



A^c = / (i'x|^c(x) 



are important because they represent constants of motion 
of the PGPE (15), and thus control the thermodynamic 
equilibrium state of the system. 



4 



3. Obtaining equilibrium properties for the C region 

To characterize the equihbrium state in the C region 
it is necessary to determine the average density, temper- 
ature and chemical potential, which in turn allow us to 
characterize the I region (see Sec. IIIB). These and other 
C region quantities can be computed by time-averaging, 
e.g, the average C region density is given by 

, Ms 

"c(x)«^ElV'c(x,t,)|', (19) 

i=i 

where {tj} is a set of Ms times (after the system has 
been allowed to relax to equilibrium) at which the field 
is sampled. We typically use 2000 samples from our sim- 
ulation to perform such averages over a time of ~ 16 
s. Another quantity of interest here is the first order 
correlation function, which we calculate directly via the 
expression 

, Ms 

G«(x,x') « — 5^ Vc(x,ti)Vc(x',t,). (20) 

i=i 

Derivatives of entropy, such as the temperature (T) 
and chemical potential (/ic) can be calculated by time 
averaging appropriate quantities constructed from the 
Hamiltonian (17) using the Rugh approach [37]. The 
detailed implementation of the Rugh formalism for the 
PGPE is rather technical and we refer the reader to Refs. 
[38, 39] for additional details of this procedure. 

A major extension to the formalism of the PGPE made 
in this work is the development of a method for extracting 
the superfluid fraction, ps, from these calculations. For 
this wo use linear response theory to relate the superfluid 
fraction to the long wavelength limit of the second order 
momentum density correlations. An extensive discussion 
of this approach, and the numerical methods used to im- 
plement it, are presented in appendix D. 



is the Hartree-Fock energy, ni is the I region density, and 
H = He + '2.h^gn\/m is the chemical potential (shifted by 
the mean-field interaction with the I region atoms) . Note 
that the average densities are constant in the uniform sys- 
tem, so Wi(k, x) has no explicit x dependence, however, 
we include this variable for generality when defining the 
associated correlation function. 

The I region density appearing in Eq. (22) is given by 

ni= / d2kW^i(k,x), (23) 

J||k||>_fC 

with corresponding atom number Ni = niL^; total num- 
ber is simply 

N = Nc + Ni. (24) 

An analytic expression for ni and simplified procedure for 
numerically calculating the first order correlation func- 
tion of the I region atoms, g[^\ can be obtained by 
taking integrals over the phase space. These results are 
discussed in appendix B. 

C. Equilibrium configurations with fixed T and TV 

Generating equilibrium classical fields with given val- 
ues of Ec and Nc is straightforward since the PGPE 
simulates a microcanonical system (see appendix A3). 
However, we wish to simulate systems with a given tem- 
perature and total number. As described in the preceding 
two sections these can only be determined after a simu- 
lation has been performed. In appendix A we outline a 
procedure for estimating values of Ec and Nc for desired 
values of N and T based on a root finding scheme using 
a Hartree-Fock-Bogoliubov analysis for the initial guess. 

IV. RESULTS 



B. Mean-field treatment of I region 

Occupation of the I region modes, A^i, accounts for 
about 25% of the total number of atoms at temperatures 
near the phase transition. We assume a time indepen- 
dent state for the I region atoms defined by a Wigner 
function [40] , allowing us to calculate quantities of inter- 
est by integrating over the above-cutoff momenta, k > K 
[41, 42]. 

Our assumed Wigner function corresponds to the self- 
consistent Hartree-Fock theory as applied in [42] . In two 
dimensions this is 

^^^'''^^= (2^)^e(^H.(k)-;l)ABT_i ' (21) 

where 

i;HF(k) = ^ + Hnc + ni), (22) 



We choose simulation parameters in analogy with the 
Paris experiment of Hadzibabic et al. [8]. This experi- 
ment used an elongated atomic cloud of approximately 
10^ *^Rb atoms, with a spatial extent (Thomas-Fermi 
lengths) of 120 /zm and 10 /zm along the two loosely 
trapped x and y directions. The tight confinement in 
the z direction was provided by an optical lattice. 

Although our simulation is for a uniform system, we 
have chosen similar parameters where possible. Our pri- 
mary simulations are for a system in a square box with 
L = 100 iim, with 4x10^ ^''Rb atoms. We also present 
results for systems with L = 50 fj.ni and L = 200 /um at 
the same density in order to better understand finite-size 
effects. All simulations arc for the case of g ~ 0.15 cor- 
responding to the experimental parameters reported in 
[8]. 

The cutoff wave vector K varied with temperature to 
ensure appropriate occupation of the highest modes (see 
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Sec. IIIAl). For the 100 system, the number of C 
region modes ranged between 559 at low temperatures to 
11338 at the highest temperature studied. 



A. Simulation of expanded interference patterns 
between two systems 

In order to make a direct comparison with the experi- 
mental results of [8] , we have generated synthetic interfer- 
ence patterns and implemented the experimental analysis 
technique. Our simulated imaging geometry is identical 
to that found in [8], with expansion occurring in the z- 
direction. The interference pattern is formed in the x-z 
plane via integration of the density along the y-direction 
( "absorption imaging" ) . 

Our algorithm for obtaining the interference pattern 
due to our classical field is very similar to that pre- 
sented in [43]. Our above cutoff thermal cloud is taken 
into account separately. We consider a pair of fields 
'ij/,^\x,y),'ip'^\x,y) from different times during the sim- 
ulation, chosen such that the fields can be considered in- 
dependent. The 3D wavefunction corresponding to each 
field is reconstructed by assuming a harmonic oscillator 
ground state in the tight-trapping direction. These two 
reconstructed fields are spatially separated by A = 3 /im, 
corresponding to the period of the optical lattice in [43] . 

Given this initial state, we neglect atomic interactions 
and only account for expansion in the tightly-trapped di- 
rection. This yields a simple analytical result for the full 
classical field ipc{x,y, z,t) at later times. The contri- 
bution of the above-cutoff atoms is included by an in- 
coherent addition of intensities. The result is integrated 
along the ^/-direction to simulate the eff'ect of absorption 
imaging with a laser beam, that is, 



nini{x,z)= dy |V'c''(x, y, z, t) |^ -|- ni(a;, y, z, t) 
Jo L 

(25) 

=^^^\x,y,z,T) + ^^^\x-A,y,z,T). (26) 

Rather than integrate the full field along the y-dircction, 
we use only a slice of length L' = 10 /xm in keeping with 
the experimental geometry of Ref . [8] . 

The interference patterns, nim(a;, z), generated this 
way contained fine spatial detail not seen in the experi- 
mental images. To make a more useful comparison to ex- 
periment it is necessary to account for the finite optical 
imaging resolution by applying a Gaussian convolution 
in the x-z plane with standard deviation 3 fim. [44]. 

In accordance with the Paris experiment, we use a 
22 ms expansion time to generate interference patterns 
for quantitative analysis (sec Sec. IV C 2). To obtain 
characteristic interference images for display in [8], the 
experiments used a shorter 11 ms expansion [44]. We ex- 
hibit examples of interference patterns at various temper- 
atures in Fig. 1, for this shorter expansion time. These 



images show a striking resemblance to the results pre- 
sented in Ref. [8]. 




20 40 
X [fim] 



20 40 
X [|im] 



FIG. 1: Synthetic interference patterns generated from the 
50 /im grid by simulation of the experimental procedure of 
Ref. [8]. (a) At low temperatures, T « O.STkt, the in- 
terference fringes arc straight, (b) Just below the transi- 
tion temperature, T w 0.95Tkt, the fringes become wavy 
due to decreased spatial phase coherence. Phase dislocations 
become common at temperatures above the transition, (c) 
T « 1.05Tkt, and (d) T I.ITkt- These "zipper patterns" 
indicate the presence of free vortices, (c) When simulation of 
the finite imaging resolution is disabled, the zipper patterns 
from the field in sub figure (d) are no longer clearly visible; the 
high-frequency details obscure the phase information without 
providing obvious additional information about the existence 
of vortex pairs. 



B. Condensate and superfluid fractions 

For a 2D Bose gas in a box we expect a nonzero con- 
densate fraction due to the finite spacing of low-energy 
modes. A central question is whether we can observe a 
distinction between the crossover due to Bose conden- 
sation and that due to BKT physics. To address this 
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question we have computed both the condensate and su- 
perfluid fractions from our dynamical simulations. 

The condensate fraction in a homogeneous system is 
easily identified as the average fractional occupation of 
the lowest momentum mode. This is directly available 
from our simulations as a time average of the k = 
mode of the classical field, 



fc = 



/N. 



(27) 



Extracting the superfluid fraction from dynamical clas- 
sical field simulations provides a more difficult challenge. 
For this we use linear response theory to relate the super- 
fluid fraction to the long wavelength limit of the second 
order momentum density correlations. Details concern- 
ing the technique are given in appendix D. 
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FIG. 2: (color online). Condensate fraction (solid dots) and 
superfluid fraction (crosses) as functions of temperature for 
the 100 grid. The transition temperature in the ther- 

modynamic limit, [14], is shown as a vertical dot-dashed 
line. The vertical dashed line shows our estimate for the tran- 
sition temperature in the finite system. The thick solid line 
is the condensate fraction for an ideal Bose gas in the grand 
canonical ensemble with the same number of atoms and pe- 
riodic spatial domain. The superfluid fraction becomes neg- 
ative in places because the extrapolation of the momentum 
correlations to k = is sensitive to statistical noise at high 
temperature (see appendix D2 for details). 

Figure 2 compares the results for the superfluid and 
condensate fractions computed on the 100 /im grid. 
These results are qualitatively similar to the results for 
the larger and smaller grids. In particular, we note that 
there is no apparent separation between temperatures at 
which the superfluid and condensate fractions fall to zero. 
Also shown in Fig. 2 is the condensate fraction for the 
ideal Bose gas confined to an identical finite-size box in 
the grand canonical ensemble. The large shift between 
ideal and computed transition temperatures indicates the 
effect of interactions in the 2D system. Because the av- 



FIG. 3: (color online). Detail of the superfluid fraction near 
the transition temperature. Solid dots represent the calcu- 
lation based on momentum correlations as described in ap- 
pendix D. Results for the largest and smallest grids are shown 
(left and right, respectively). The data for the 100 fim grid 
is omitted for clarity, but lies between the curves shown as 
expected. Open circles represent the calculation of the super- 
fluid fraction from the associated fitted values for the decay 
coefficient a, via Eq. (4). The open circles terminate where 
the power law fitting procedure fails. 



erage system density is uniform, this large shift is to due 
to critical fluctuations (also see [34]). 

In our calculations we identify the transition tempera- 
ture, Tkt, as where the superfluid fraction falls off most 
rapidly (i.e., the location of steepest slope on the ps ver- 
sus T graph; see Fig. 2). As the system size increases, 
this transition temperature moves toward the value for an 
infinite-sized system, [14]. This effect is illustrated 
by the behavior of the superfluid fraction in Fig. 3. 



C. First order correlations — algebraic decay 

Algebraic decay of the first order correlations, as de- 
scribed by Eq. (3), is a characteristic feature of the BKT 
phase. Above the BKT transition, the first order corre- 
lations should revert to the exponential decay expected 
in a disordered phase. 

The normalized first order correlation function, 5(1) is 
defined by 



gW(x,x') = 



G(i)(x,xO 
■\/n(x)ri(x') ' 



(28) 



where G^^^{x,x.') = (■!/;'l'(x)'(/'(x')) is the unnormalized 
first order correlation function [40]. 
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1. Direct calculation of g^^^ 

In the PGPE formalism the C and I contributions to 
the correlation function are additive [41], that is, 



G(i) (x, x') = Gg^ (x, x') + G['' (x, x' 



(29) 



where Gq^ and G^' are defined in Eqs. (20) and (B8), 

respectively. It is interesting to note that Gq'' and Gj^'* 
individually display an oscillatory decay behavior — orig- 
inating from the cutoff — an effect which correctly can- 
cels when the two are added together. 

Having calculated g^^\ we obtain the coefficient a by 
fitting the algebraic decay law, Eq. (3), using nonlinear 
least squares; sample fits are shown in Fig. 4. The fit is 
conducted over the region between 10 and 40 de Broglie 
wavelengths. The short length scale cutoff is to avoid 
the contribution of the non-universal normal atoms, for 
which the thermal de Broglie wavelength sets the appro- 
priate decay length. The long distance cutoff is chosen 
to be small compared to the length scale L, to avoid the 
effect of periodic boundary conditions on the long range 
correlations. 

The quality of the fitting procedure, and the break- 
down of expression (3) at the BKT transition can be ob- 
served by adding an additional degree of freedom to the 
fitting function. In particular, at each temperature we 
fit the quadratic hi(g^^^ ) = A — a hi{x) + S In^ (x) and ex- 
tract the parameter S {a ~ a is discarded). The abrupt 
failure of the fits can be observed in the inset of Fig. 5 
as a sudden increase in the value of |^(r)| — an effect 
which is in excellent agreement with the value of Tkt as 
estimated from the superfluid fraction. 
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FIG. 4: (color online). Sample fits to the algebraic decay 
of 5*1) at various temperatures, ranging from below to above 
the transition. High temperatures correspond to curves at 
the bottom of the figure which have rapid falloff of g'^' with 
distance. Fits are shown on a log-log scale in the inset to 
emphasize the failure of a power law in describing the behavior 
of g^-^^ at high temperature. 
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0.1 



50 100 150 200 

T[nK] 

FIG. 5: (color online). Comparison of two methods for deter- 
mining the algebraic decay coefficient a{T) for the first order 
correlation function (;'^'(x,x'). The line with circle markers 
represents direct fits to g'^'. These fits fail at the transition 
temperature as shown by the sharply diverging value of \S{T) 
in the inset. The filled points represent the values a' (T) ob- 
tained from a simulation of the experimental analysis proce- 
dure of [8], described in Sec. IV C 2. Horizontal dotted lines 
at 0.25 and 0.5 correspond to the expected values of a' just 
below and above the transition, respectively [8]. The vertical 
line is the BKT transition temperature, as estimated from the 
superfluid fraction calculated in Sec. IV B. 



2. Calculation of g^^^ via interference patterns 

So far a direct probe of the in situ spatial correlations 
has not been possible, although important progress has 
been made by the NIST group [11]. In the experiments of 
Hadzibabic et al. [8] a scheme proposed by Polkovnikov 
et al. [45] was used to infer these correlations from the 
"waviness" of interference patterns produced by pair of 
quasi-2D systems (see Sec. IV A). In this section we simu- 
late the experimental data analysis method, and compare 
inferred predictions for the correlation function against 
those we can directly calculate. This allows us to char- 
acterize the errors associated with this technique arising 
from finite size effects and finite expansion time. 

To make this analysis we follow the procedure outlined 
in [8] . We fit our numerically generated interference pat- 
terns (see Sec. IV A) to the function 



F(x,z) = G{z) 



1 -f c{x) COS 



27rz 



-e{x] 



(30) 



where G{z) is a Gaussian envelope in the z-direction, c{x) 
is the interference fringe contrast, D is the fringe spacing 
and 0{x) is the phase of the interference pattern in the 
z-direction. 

Defining the function 

C{L^) = ^ c{x)e''^^Ux, (31) 



ix/2 
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the nature of spatial correlations is then revealed by 
the manner in which <^|C(La;)p) decays with Lj.. In 
particular, we identify the parameter a', defined by 
{\C[L^)\^) cx [45]. For an infinite 2D system in the 

superfluid regime (T < T^) a' = a (i.e. a' corresponds 
to the algebraic decay of correlations). For T > T^, 
where correlations decay exponentially, a' is equal to 0.5. 

Fitting (|C(La;)p) to the algebraic decay law AL'"^""' 
we can determine a'. A comparison between a' inferred 
from the interference pattern and a obtained directly 
from g(^) is shown in Fig. 5. Both methods give broadly 
consistent predictions for a when T < TkT) however our 
results show that there is a clear quantitative difference 
between the two schemes, and that a' underestimates the 
coefficient of algebraic decay in the system (i.e. using a' 
in Eq. (4) would overestimate the superfluid density). 
Near and above transition temperature, where the fits to 
g^^^ fail, we observe that a' converges toward 0.5. The 
agreement between a and a' in the low temperature re- 
gion improves as the size of the grid is increased. 



D. Vortices and pairing 

The simplest description of the BKT transition is that 
it occurs as a result of vortex pair unbinding: At T < Tkt 
vortices only exist in pairs of opposite circulation, which 
unbind at the transition point to produce free vortices 
that destroy the superfluidity of the system. However, to 
date there are no direct experimental observations of this 
scenario, and theoretical studies of 2D Bose gases have 
been limited to qualitative inspection of the vortex dis- 
tributions. In the c-field approach vortices and their dy- 
namics are clearly revealed, unlike other ensemble-based 
simulation techniques where the vortices are obscured by 
averaging. This gives us a unique opportunity to inves- 
tigate the role of vortices and pairing in a dilute Bose 
gas. 

We detect vortices in the c-field microstates by ana- 
lyzing the phase profile of the instantaneous field (see 
appendix C). An example of a phase profile of a field 
for T < Tkt is shown in Fig. 6(a). The vortex locations 
reveal a pairing character, that is, the close proximity of 
pairs of positive (clockwise) and negative (counterclock- 
wise) vortices relative to the average vortex separation. 
An important qualitative feature of our observed vortex 
distributions is that at high temperatures, pairing does 
not disappear from the system entirely. Indeed, most vor- 
tices at high temperature could be considered paired or 
grouped in some manner, as shown in Fig. 6(b). Perhaps 
this is not surprising, since positive and negative vor- 
tices have a logarithmic attraction, and we observe them 
to create and annihilate readily in the c-field dynamics. 
However, this does indicate that the use of pairing to lo- 
cate the transition may be ambiguous, and we examine 
this aspect further below. 

It is also of interest to measure the number of vortices, 
Ny, present in the system as a function of temperature 
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FIG. 6: (color online). Pliase profile of a c-field with vortices 
indicated. Vortices with clockwise (white -|-) and anticlock- 
wise (black o) circulation. The phase of the classical field 
is indicated by shading the background between dark blue 
(phase 0) and light yellow (phase 2tt) . (a) Distinctive pairing 
below the transition at T = 207nK ^ 0.93rKT (b) A "vortex 
plasma" above the transition at T = 238nK ~ I.OTTkt. 



(see Fig. 7). At the lowest temperatures the system is 
in an ordered state, and the energetic cost of having a 
vortex is prohibitive. As the temperature increases there 
is a rapid growth of vortex population leading up to the 
transition point followed by linear growth above Tkt- 



1. Radial vortex density 

The most obvious way to characterize vortex pairing 
is by defining a pair distribution function for vortices of 
opposite sign. Adopting the notation of [26], this is 

Gi%(r) = (p„,+ (0)p,„,_(r)), (32) 
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FIG. 7: Total number of vortices (dots) and number of un- 
paired vortices (circles) as a function of temperature near the 
transition. While at the transition temperature is already 
very high, Nu becomes nonzero only close to the transition, 
providing clear evidence of vortex unbinding at work. The 
inset shows the variation in the total number over the full 
temperature range of the simulations. Above the transition 
temperature the growth in the number of vortices becomes 
linear with temperature. 



where ^ is the vortex density function which consists 
of a sum of delta-spikes, 



E 



S{r-vt) 



for positive vortices at positions {fi^}- We use the anal- 
ogous definition for Pv,-- The associated dimensionless 
two- vortex correlation function is 



(r) 



9v,± 



(r) 



(p,,+ (0))(p,,_(r)) 



(33) 



(2) 

The angular average of g„ 4 can be calculated directly 
from the detected vortex positions using a binning pro- 
cedure on the pairwise distances ||r^ — r~ || , and is shown 
in Fig. 8. 

These results quantify the effect discussed earlier: Pos- 
itive and negative vortices show a pairing correlation 
which does not disappear above Tkt- The characteristic 
size of this correlation, given by twice the width of the 
peak feature in Fig. 8, is Icoi ^ 3/im (taking full width 
half maximum) . 

The shape of our pairing peak is qualitatively similar 
to that described in [26]. However, in contrast to their 
results the width does not appear to change apprecia- 
bly with temperature. Additional simulations show that 
increasing the interaction strength causes the peak to be- 
come squarer and wider. It is clear that while the pair 

(2) 

size and strength revealed in 5^ ±{r) does not change 
appreciably as the transition is crossed, the amount of 
pairing relative to the background uncorrelated vortices 



changes considerably. This background of uncorrelated 

(2) 

vortices is given by the horizontal plateau ^.{r) ^ 1 at 
large r as shown in the inset. 




FIG. 8: (color online). Angular average of the two- vortex pair 
distribution functions for vortices of opposite sign. Three 
temperatures centered about the transition are shown: dot 
markers T = 194nK ^ O.gTKT, fc = 0.34; circle markers 
T = 217nK ^ I.OITkt, fc = 0.076; cross markers T = 236nK 
~ I.ITkt, fc ~ 0.006. The vertical dotted line shows the 
value of the healing length at T = 0. The main plot shows 
gl^l- normalized by the positive vortex density; comparable 
magnitudes for the peaks near r = show that vortex pairing 
remains important over the range of temperatures studied, 
not only below the transition. The inset shows g^l- in the 
natural dimensionless units for which gl^±{r) — >■ 1 as r — > 00. 



2. Revealing unpaired vortices with coarse- graining 
(2) 

The function GJ, 4 (r) clearly indicates the existence 
of vortex pairing in the system. However, it does not 
provide a convenient way to locate the transition tem- 
perature, since a large amount of pairing exists both 
below and above the transition: The expected number 
of neighbors for any given vortex — roughly, the area 

(2) 

of the pairing peak of {ny^+)Gl !^{r) shown in Fig. 8 
— does not change dramatically across the transition. 
(ni,^+) = (rii,)/2 is the expected density of positive vor- 
tices. 

We desire a quantitative observation of vortex unbind- 
ing at the transition and have therefore investigated sev- 
eral measures of vortex pairing ^. However, measures 
based directly on the full set of vortex positions seem to 
suffer from the proliferation of vortices at high tempera- 



^ For example, the Hausdorff distance (see, e.g., [46]) between the 
set {r^} of positive vortices and the set {r~ } of negative vortices. 
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ture — an effect which tends to wash out clear signs of 
vortex unbinding. With this in mind, we have developed 
a procedure for measuring the number of unpaired vor- 
tices in our simulations, starting from the classical field 
rather than the full set of vortex positions. 

The basis of our approach for detecting unpairing is 
to coarse-grain the classical field by convolution with a 
Gaussian filter of spatial width (standard deviation) aj. 
This removes all vortex pairs on length scales smaller 
than than af. Figure 9 shows the count of remaining 
vortices as a function of filter width, along with some ex- 
amples of coarse-grained fields. For (Jf > Icor, the number 
of remaining vortices levels off and only decreases slowly 
with increasing af. Ultimately the number of remaining 
vortices goes to zero as af ^ L. 

Setting the filter width to be larger than the charac- 
teristic pairing distance, Icor, yields a coarse-grained field 
from which the pairs have been removed, but unpaired 
vortices remain. In our simulations we have /cor ~ 3 /im; 
we take the vortices which remain after coarse-graining 
with a Gaussian of standard deviation crj = 5 /im to 
give an estimate of the number of unpaired vortices, Nu- 
Figure 7 shows that Nu becomes nonzero only near the 
transition, in contrast to Ny which is nonzero well below 
Tkt- The sharp increase in Nu at Tkt is a quantitative 
demonstration of vortex unbinding at work. 
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FIG. 9: (color online). The coarse-graining procedure: num- 
ber of vortices as a function of filter width for a tempera- 
ture near the transition. The smooth curve is an average 
over many realizations of the field, whereas the stepped curve 
shows typical behavior of the number for a single field. In- 
sets show the coarse-grained fields for various filter widths; 
the transformation removes vortex- ant i vortex pairs which are 
separated by approximately less than the standard deviation 
of the filter. In this example Nu = 4 unpaired vortices remain 
at o"/ = 5 /im. 

In the experiment of Ref. [8] , the fraction of interfer- 
ence patterns with dislocations (e.g., see Figs. 1(c) and 
(d)) was measured. While isolated vortices are clearly 
identified by interference pattern dislocations, a lack of 



spatial resolution in experiments means that this type 
of detection method obscures the observation of tightly 
bound vortex pairs. The experimental resolution of 3 /xm 
is broadly consistent with the scale of the coarse-graining 
filter (i.e., af = 5 /xm). With this in mind, we introduce 
the quantity pu (T) , defined as the probability of observ- 
ing an unpaired vortex in a 50 x 50 /im control volume 
at a given temperature ^. For the 50 /im grid we have 
simply Pu{T) = Pt{Nu > 1). 

Computing Pu{T) from our simulations yields the re- 
sults shown in Fig. 10. Our results show a dramatic jump 
in Pu at a temperature that is consistent with the tran- 
sition temperature Tkt determined from the superfiuid 
fraction calculation presented in Sec. IV B. 
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FIG. 10: (color online) . Comparison of vortex unpairing mea- 
sures. The dots are our pairing measure based on coarse- 
graining the field. Circles represent the pairing as determined 
by the number of dislocations in the simulated interference 
patterns. This was the same method used in the experimental 
analysis of [8] and coincides remarkably well with our coarse- 
graining based measure. Both curves are consistent with the 
vertical line showing the transition temperature Tkt as deter- 
mined from the superfiuid fraction calculation in Sec. IV B. 
The inset shows the calculated coarse-grained pairing mea- 
sure for all three grid sizes, along with vertical lines showing 
the estimates for Tkt derived from the superfiuid fraction 
calculations. 

From the definition, we expect that p„ should be close 
to the experimentally measured frequency of dislocations. 
To demonstrate this relationship, we have simulated in- 
terference patterns (as described in Sec. IV A) and de- 
tected dislocations using the experimental procedure of 
Ref. [8]: A phase gradient dO/dx was considered to mark 
a dislocation whenever \d9/dx\ > 7r/4rad//im. From this 
we can compute the probability of detecting at least one 
dislocation as a function of temperature. As shown in 



* We choose a fixed control volume with L = 50 /im in order to 
compare results between simulations with different grid sizes. 
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Fig. 10, the results of this procedure compare very fa- 
vorably with our measure of pairing based on We 
note that inhomogeneous effects in experiments proba- 
bly broaden the jump in p„ appreciably compared to our 
homogeneous results. 



V. CONCLUSION 

In this paper wc have used c-field simulations of a 
finite-sized homogeneous system in order to investigate 
the physics of the 2D Bose gas in a regime corresponding 
to current experiments. We have directly computed the 
condensate and superfluid fractions as a function of tem- 
perature, and made comparisons to the superfluid frac- 
tion inferred by the first order correlation function, and 
using the interference scheme used in experiments. Our 
results for these quantities provide a quantitative test of 
the interference scheme for a finite system. 

A beautiful possibility is the direct experimental ob- 
servation of vortex- antivortcx pairs, their distribution in 
the system, and hence a quantitative measurement of 
their unbinding at the BKT transition. We have calcu- 
lated the vortex correlation function across the transition 
and provided a coarse-graining scheme for distinguishing 
unpaired vortices. These results suggest that the dislo- 
cations observed in experiments, due to limited optical 
resolution, provide an accurate measure of the unpaired 
vortex population and accordingly are a strong indicator 
of the BKT transition. 

We briefly discuss the effect that harmonic conflnement 
(present in experiments) would have on our predictions. 
The spatial inhomogeneity will cause the superfluid tran- 
sition to be gradual, occurring first at the trap center 
where the density is highest, in contrast to our results 
where the transition occurs in the bulk. So far the super- 
fluid fraction for the trapped system has been determined 
by using the universality result for the critical density in 
the homogeneous gas [15] in combination with the local 
density approximation [21, 25]. It would be interesting to 
be able to compute the superfluid fraction independently 
as we have done here; however it is not clear how to do 
so. 

Bisset et al. [25] used an extension of the c-field method 
for the trapped 2D gas to examine g^^^ and found similar 
results for the onset of algebraic decay of correlations at 
the transition. Their analysis was restricted to the small 
region near the trap center where the density is approxi- 
mately constant; we expect the results of our vortex cor- 
relation function and the coarse-graining scheme should 
similarly be applicable to the trapped system in the cen- 
tral region. Except in very weak traps, the size of this 
region is relatively small and will likely prove challenging 
to measure experimentally. 

Our results for the homogeneous gas emphasize the 
clarity with which ab initio theoretical methods can cal- 
culate quantities directly observable in experiments, such 
as interference patterns. This should allow direct com- 



parisons with experiments, providing stringent tests of 
many-body theory. 
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Appendix A: Simulation using the PGPE 

Here we outline our procedure for determining the 
properties of the C region and the steps used to cre- 
ate initial states for the PGPE solver. The C region 
itself is characterized by the cutoff momentum K, while 
the initial states are characterized by the energy Ec and 
number A^c- We want to obtain values of these three 
properties which are consistent with a specified temper- 
ature T and total number of atoms N. 



1. HcU-tree-Fock-Bogoliubov analysis 

To generate an initial estimate of the C region parame- 
ters we solve the self-consistent Hartree-Fock-Bogoliubov 
(HFB) equations in the so-called Popov approximation 
[47] to find an approximate thermal state for the system 
at a temperatm'e T. The resulting state is a Bose Ein- 
stein distribution of quasiparticles interacting only via 
the mean-field, expressed in terms of the quasiparticle 
amplitudes w,k and Wk- 

Occupations for the C region field may be computed 
directly from the quasiparticle occupations via 

ni,= {ul + vl)NB{Ei,)+vl, (Al) 

where Nb is the Bose Einstein distribution and is 
the quasiparticle energy which is obtained by solving the 
Bogoliubov-de Gennes equations self consistently [47]. 
This allows us to compute the cutoff as the maximum 
value of ||k|| consistent with sufficient modal occupation: 

/s: = max{||k||: rik > ncut}- (A2) 

We choose ricut = 5 for the sufficient occupation condi- 
tion on the C region modes. 

The number of atoms below the cutoff may be com- 
puted directly from the sum of the condensate number 
A^o and the number of C region excited state atoms, Nic' 

Nc=No+Nic, where A^ic = ^ nk- (A3) 

kec\{o} 
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For the total energy below the cutoff, we use the ex- 
pression 



^ E^[Nb{E^)-vI] (A4) 
kec\{o} 



where A = g{NQ + 2Nic)- Rearranging, this is 



kGC\{0} 

(A5) 



The expression (A5) differs from Eq. (22) of [47] as we 
have retained the zcroth order (constant) terms that are 
required to match the energy scale of the HFB analysis to 
the zero point of energy in the classical field simulations. 



2. Initial conditions for fixed total number 

A simple comparison between simulations at varying 
temperatures can only be carried out if the total number 
of atoms is fixed. This presents a problem in our sim- 
ulations: although the number of atoms and energy of 
the C region can be directly specified (see Sec. A3), we 
may only determine the total number after performing a 
simulation. This is because the number of atoms in the I 
region depends on the temperature and chemical poten- 
tial which are calculated by ergodic averaging of the C 
region simulations. 

Formally, this may be stated as a root-finding problem: 
solve 



N{Nc,Ec) 



(A6) 



with initial guess provided by the solution to the HFB 
analysis in Sec. Al. Although both A^c and Ec affect 
the total number N, we choose to fix A^c to the initial 
guess and to vary Ec until the desired total number is 
found. 

We note that evaluating the function N{Nc, Ec) is 
very computationally expensive and difficult to fully au- 
tomate since it involves a simulation and several steps of 
analysis. For this reason we use a nonstandard root find- 
ing procedure: For the first iteration we simulate three 
energies about the initial guess Ec such that the results 
crudely span A^tot; these three simulations can be per- 
formed in parallel which significantly reduces the time to 
a solution. A second guess was obtained by quadratic 
fitting of Ec as a function of N which gives N accurate 
to within about 5% of A'tot- An addition iteration using 
the same interpolation method takes N to within 0.3%, 
which we consider sufficient. 

We note that changing Ec during the root finding pro- 
cedure means we have no direct control over the final 
temperature of each specific simulation. In our case this 



is not a problem because we only require a range of tem- 
peratures spanning the transition. In principle one could 
solve for a given temperature by allowing A^c to vary in 
addition to Ec- 



3. Initial conditions for given Ec and Nc 

We compute initial conditions for the C region field in 
a similar way to [48] . Using the representation for the C 
region given by Eq. (10), the task is to choose appropriate 
values for the {c„}. As a first approximation, choose the 
smallest value for a momentum cutoff K' such that the 
field with coefficients 



Cn = 




for < ||k|| < K', 
for Ikl > K\ 



(A7) 



has energy greater than Ec- Here A is chosen so that 
the field has normalization corresponding to A'c atoms, 
and 0n is a randomly chosen phase which is fixed for each 
mode at the start of the procedure. The random phases 
allow us to generate many unique random initial states 
at the same energy. 

By definition, the field defined by (AT) has energy 
slightly above the desired energy. This problem is solved 
by mixing it with the lowest energy state: 




for n = 0, 

elsewhere, 



(A8) 



using a root finding procedure to converge on the desired 
energy Ec- The scheme generates random realizations 
of a non-equilibrium field with given Ec and Nc which 
are then simulated to equilibrium before using ergodic 
averaging for computing statistics. 



Appendix B: I region integrals 

Our assumed self-consistent Wigner function 
(Sec. IIIB) for the I region atoms takes a particu- 
larly simple form in the homogeneous case: 



1 



(27r)2 g(R2k2/2m-|-2B2g„c/TO-Mc)/fcBT _ l' 

(Bl) 

The above-cutoff density may then be found by direct 
integration: 



ni(x) = / d^k Wi(k, 

J\M>K 



(B2) 



_ ^-(h'^K'^ /2m+2h^gnc/m-nc)/kBT 

(B3) 

with A the thermal de Broglie wavelength. 
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In a similar way, the assumed Wigner function allows 
any desired physical quantity to be estimated via a suit- 
able integral. A particiilar quantity of interest in the 
current work is the first order correlation function, which 
can be obtained from the Wigner function as [40] 

Gi^^(x,x') = / rf^k e-'''-("-"') Wi(k, 

J\\k\\>K 

(B4) 

This integral is of the general form 

for constants A and B. Noting that Ii depends only on 
the length, r of ||r||, and transforming k to polar coordi- 
nates {k,6), we have 

/oo /•2tt 
^ ^ %A.^+B_i I rf^e-™^ (B6) 

= ^ %A«^+B_i 2[r(i)]Vo(m), (B7) 

(see [49, p902] for the Bessel function identity). 

Thus we obtain g[^\x,x') in terms of a one dimen- 
sional integral which may be performed numerically: 




kJo(«;||x — x'll) 

.{H^K^/2m+2h^gnc/m-iJ,c)/kBT _ l' 



(B8) 

Appendix C: Vortex detection 

The defining feature of a "charge-m" vortex is that the 
phase 6 of the complex field '4'{x) = |'!/;(x)|e'^*^''-* changes 
continuously from to 2rmr around any closed curve 
which circles the vortex core. We express our field tp 
on a discrete grid in position space; the aim of vortex de- 
tection is then to determine which grid plaquettes (that 
is, sets of four adjacent grid points) contain vortex cores. 

To obtain the phase winding about a plaquette, first 
consider the phase at two neighboring grid points A and 
B. We are interested in the unwrapped phase difference 
A^AB between the grid points; unwrapping ensures that 
the phase is continuous between A and B. (In the discrete 
setting such continuity is poorly defined; the best we can 
do is to correct for the possibility of 2tt phase jumps by 
adding or subtracting factors of 2tt so that |A^ab| < tt.) 
The unwrapped phase differences around a grid plaquette 
tell us a total phase change ^wrap = X^i = 2rmr 

where m £ Z is the winding number or "topological 
charge" . 

Due to the necessity of unwrapping the phase, a four- 
point grid plaquette cannot unambiguously support vor- 
tices with charge larger than one. Luckily, such vortices 
are energetically unfavorable in 2D Bose gases [50, p. 83] 
so we need only concern ourselves with detecting vortices 



with winding number ±1 in this work. The positions ob- 
tained from a given run of our vortex detection algorithm 
are the labeled {r^*"} and {r^} for winding numbers -|-1 
and —1, respectively. 



Appendix D: Superfluid fraction 

One of the important characteristics of the BKT tran- 
sition is the presence of superfluidity, even in the absence 
of conventional long range order. In the following we de- 
scribe a method to calculate the superfluid fraction from 
our classical field description. The method is attractive 
because it makes use of momentum correlations which 
may be extracted directly from our equilibrium simula- 
tions without any need to introduce additional boundary 
conditions or moving defects. 

1. Superfluid fraction via momentum density 
correlations 

Our derivation is based on the procedure presented in 
[51, p.214], (see also [52] and [50, p.96]). The central 
idea is to establish a relationship between i) the auto- 
correlations of the momentum density in the simulated 
ensemble and ii) the linear response of the fluid to slowly 
moving boundaries; (i) is a quantity we can calculate, 
while (ii) is related to the basic properties of a superfluid 
via a simple thought experiment. 

To connect the macroscopic, phenomenological de- 
scription of supcrfliiidity with our microscopic theory, we 
make use of the thought experiment shown schematically 
in Fig. 11(b): Consider an infinitely long box, B contain- 
ing superfluid, and accelerate the box along its long axis 
until it reaches a small velocity u. Due to viscous inter- 
actions with the walls, such a box filled with a normal 
fluid should have a momentum density at equilibrium of 
(p)u = nu. The notation denotes an expectation 
value in the ensemble with walls moving with velocity u. 

Because superfluids are nonviscous, the observed value 
for the momentum density in a superfluid is less than the 
value nu expected for a classical fliiid. In Landau's two- 
fluid model we attribute the observed momentum density, 
p„u, to the "normal fraction" where p„ is the normal 
fluid density. The superfluid fraction remains stationary 
in the lab frame, even at equilibrium and makes up the 
remaining mass with density ps = n — pn- 

In order to apply the usual procedures of statistical 
mechanics to the thought experiment, we consider two 
frames: the "lab frame" in which the walls move with ve- 
locity u in the x-direction and the "wall frame" in which 
the walls are at rest. 

Assuming that the fluid is in thermal eqiiilibrium 
with the walls, the density matrix in the grand canon- 
ical ensemble is given by the usual expression p = 
^-f}{H^-^N) I rj,^ l^^-fi{H^-i,N)-^ where is the Hamil- 

tonian of the system in the wall frame and = l/ksT. 
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FIG. 11: Thought experiment used in deriving the superfluid 
density. The walls move with velocity u in the x direction. 
To begin with, we imagine that the superfluid sits in a box of 
dimensions x Ly as shown in a). We later take the limit 
as the box walls recede to infinity to get the thermodynamic 
limit d). The order of the limits is critically important: the 
path b) leads to supcrfiow while the path c) results in the 
entire fluid moving along with the walls. 

A Galilean transformation relates to the Hamilto- 
nian in the lab frame, = iJ — u • P + |Mw^, where 
P = Jg p(x)(i^r is the total momentum, M = niN is the 
total mass and p(x) is the momentum density operator 
at point X. The expectation value for the momentum 
density in the presence of moving walls is then given by 
the expression 

(p(x))„ = Tr[pp(x)], (Dl) 
= — ^ . (D2) 

Expanding this expression to first order in u yields 

(p(x))„ = (p(x))+/3((p(x)P.u)-(p(x)) (P.u)), (D3) 

where all the expectation values on the right hand side 
are now taken in the equilibrium ensemble with the walls 
at rest. Since (p(x)) = in our equilibrium ensemble, 
this simplifies to 

(p(x))„ = /3(p(x)P) • u, (D4) 

= 13 [ dP^ (p(x)p(x')) • u, (D5) 

JB 

where p(x)p(x') is a dyad [i.e., a rank-two tensor; the 
outer product of p(x) and p(x')]. 

To make further progress, wc wish to take the limit 
as the system gets very large (write "i? — >■ oo"). To 
this end, we first consider some properties of the corre- 
lation functions in the infinite system. The infinite sys- 
tem is homogeneous, which implies that (p(x)p(x'))^ = 



(p(x + r)p(x' -I- r))^ for any r, where {■) ^ indicates an 
average in the infinite system. As a consequence, we may 
express the correlations — in the infinite system — in 
terms of the Fourier transform in the relative coordinate 
x' — x: 

(p(x)p(x'))^ = (p(0)p(x'-x))^, (D6) 
= (2^ /^i'ke^'-(-'--)x(k), (D7) 

where all the important features of the correlations are 
now captured by the tensor 

x(k)=y'd2re-^'^-(p(0)p(r))^. (D8) 

Because of the isotropy of the fluid in the infinite system, 
x(k) obeys the transformation law x(Ok) = 0~^x(k)0 
for any 2x2 orthogonal matrix O. This implies that 
X may be decomposed into the sum of longitudinal and 
transverse parts: 

X(k) =kkxi(fc) + (/-kk)xt(fc) (D9) 

where k = k/fc,A;=||k|| and / is the identity. The trans- 
verse and longitudinal functions Xt and xi are scalars 
which depend only on the length k. 

We now return our attention to the finite system. If 
the finite box B is large then the momentum correlations 
in the bulk will be very similar to the values for the infi- 
nite system. Therefore, when x and x' are far from the 
boundaries, we may approximate 

(p(x)p(x')) « (P(x)p(x')),, (DIO) 

= (2^ |d'ke*-('''--)x(k) (Dll) 

which in combination with Eq. (D5) yields 

(p(x))„«/3 j/^J^ I rf^ke*-(-'--)x(k)-u 

(D12) 

= ^ y d^k As (k)e* ''x(k) • u. (D13) 

Here we have defined the nascent delta function 
AB(k) :— Jg (Px' e^^'^ which has the property 

AB(k) (5(k) as S -j> oo. 

We are now in a position to carry out the limiting 
procedure to increase the box size to infinity. How- 
ever, care must be taken because the simple expression 
lirns^oo (p(x))^ is not well defined without further qual- 
ification of the limiting process B ^ oo. 

To resolve this subtlety we must insert a final vital 
piece of physical reasoning. Let us assume for simplicity 
that u is directed along the x-direction, and the box B is 
aligned with the x and y axes with dimensions L^xLy. 
As shown in Fig. 11, there are two possibilities for taking 
the limits, representing different physical situations. 
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On the one hand [Fig. 11(b)], we may take the Umit 
Lx oo first, which gives us an infinitely long channel 
in which superfluid can remain stationary while only the 
normal fraction moves with the walls in the x-direction. 
We have 



p„u = lim lim (p(x))^, 



(DM) 



lim lim /3 /(i2kAs(k)e*-^x(k)-u (D15) 

(D16) 



= j3 lim lim x(k) • u 



where we use the fact that AB(k) can be decomposed 
into the product /S.L^{kx)^Ly{ky) with Ai(fc) 5{k) 
as L — > oo. Employing the decomposition of x given in 
Eq. (D9) allows the density of the normal fraction to be 
related to the transverse component of x evaluated at 
zero: 



' lim Xt{k) 



/3Xt(0). 



(D17) 



On the other hand [Fig. 11(c)] we may take the limit 
Ly ^ oo first, resulting in an infinitely long channel — 
with velocity perpendicular to the walls in which the 
entire body of the fluid must move regardless of the su- 
perfluidity. In a similar way to the previous paragraph, 
nu = /31imfe^^o limfe^^o x(k) • u, and making use of the 
decomposition in Eq. (D9) , the total density is related to 
the longitudinal component of the correlations: 



n = /3 lim Xi(fc) = /3x;(0)- 

k->-0 



(D18) 



With these expressions, the normal fraction /„ may 
finally be expressed directly as 



fn = Pn/n ■■ 



lim Xt(A:)/ linixi(^) 

fe->-0 fe->0 



(D19) 



while the superfluid fraction is fa ~ 1^ fn- Thus, we have 
expressed the superfluid and normal fractions in terms of 
a correlation function which can be directly computed 
from our simulation results. 



2. Numerical procedure 

To determine the superfluid fraction for our system, 
we need to estimate the tensor of momentum density 
correlations x from the simulation results. For our flnite 



system constrained to a periodic simulation box, we may 
compute the momentum correlations only at discrete grid 
points. The discrete analogue of Eq. (D8) leads to the 
expression 



X(k) oc (pkP-k) 



(D20) 



where the constant of proportionality is not important 
to the final result, and pk are the discrete Fourier coeffi- 
cients of p(x) over our simulation box. 

The momentum density operator is given by 



p(x) = I [(V^t(x))^(x) 



V'1'(x)VV'(x) (D21) 



which may be derived by considering the continuity equa- 
tion for the number density, ('0^(x)^/;(x)). For a given 
classical field Eq. (10), the Fourier coefficients of p may 
be written as 



Pk 



^(2k' + k)c^,Ck+k' 



(D22) 



where As is the area of the system. Computing a value 
for all Pk at each time step, we then evaluate x(k) via 
the usual ergodic averaging procedure using Eq. (D20). 

Having evaluated x(k), we are left with performing 
the decomposition into longitudinal and transverse parts. 
For this, simply note that Eq. (D9) implies xi{k) = k • 
x(k)-k, and Xt{k) = w-x(k)-w, where w is a unit vector 
perpendicular to k. 

Values for Xt and xi may be collected for all angles as 
a function of k, and a fitting procedure used to perform 
the extrapolation fc 0; this procedure is illustrated 
in Fig. 12. At low temperatures, the extrapolation is 
quite reliable, but becomcis more difficult near the tran- 
sition where sampling noise increases and Xt{k) changes 
rapidly near fc = 0. Without a known functional form, 
we settled for a quadratic weighted least squares fit of 
In(xt) and In(xi) versus k. A weighting of 1/k was used 
to counteract the fact that the density of samples of x 
vs fc scales proportionally with fc due to the square grid 
on which x(k) is evaluated. The logarithm was used to 
improve the fits of Xt very near the transition where it 
varies non-quadratically near fc = 0. The fitting proce- 
dure and extrapolation to fc = generally produces rea- 
sonable results, but is somewhat sensitive to numerical 
noise. For this reason, the computed superfluid fraction 
at high temperatures is not exactly zero (see Fig. 2). 
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